*
* $Id$
*
* $Log: gmicap.F,v $
* Revision 1.1.1.1  2002/06/16 15:18:43  hristov
* Separate distribution  of Geant3
*
* Revision 1.1.1.1  1999/05/18 15:55:21  fca
* AliRoot sources
*
* Revision 1.1.1.1  1995/10/24 10:21:52  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.47  by  S.Giani
*-- Author :
      SUBROUTINE GMICAP
C
#include "geant321/gctrak.inc"
#include "geant321/gcmate.inc"
#include "geant321/gcking.inc"
C MICAP commons
#include "geant321/mmicap.inc"
#include "geant321/minput.inc"
#include "geant321/mconst.inc"
      COMMON/MNUTRN/NAME,NAMEX,E,EOLD,NMED,MEDOLD,NREG,U,V,W,
     + UOLD,VOLD,WOLD,X,Y,ZZ,XOLD,YOLD,ZOLD,WATE,OLDWT,WTBC,
     + BLZNT,BLZON,AGE,OLDAGE,INEU,ENE(MAXNEU)
      INTEGER BLZNT
#include "geant321/mapoll.inc"
#include "geant321/mpoint.inc"
#include "geant321/mrecoi.inc"
#include "geant321/mmass.inc"
#include "geant321/mpstor.inc"
#include "geant321/cmagic.inc"
#include "geant321/mcreco.inc"
C
C convert Z,A of recoil to CALOR particle code
C only p = 0, D = 7, T = 8, He3 = 9, alpha=10
* n = 13, p = 14, D = 45, T = 46, He3 = 49, alpha = 47
      DIMENSION NGPART(4,0:2)
      DATA ((NGPART(I,J),I=1,4),J=0,2)/13 ,-1 ,-1 ,
*                                      n
     +                               -1 , 14 , 45 ,
*                                         p    D
     +                               -1 , 46 , 49 ,
*                                         T    He3
     +                               -1 ,-1  , 47/
*                                              alpha
*      SAVE
C first check, if ZEBRA still in order
      IF(LD(LMAG1).NE.NMAGIC.OR.LD(LMAG2).NE.NMAGIC) THEN
         WRITE(6,*) ' CALOR: ZEBRA banks screwed up --> STOP'
         WRITE(IOUT,'('' MICAP: Magic number '',I12,'' not found: '',  '
     +   //'      2I12)') NMAGIC,LD(LMAG1),LD(LMAG2)
         STOP
      ENDIF
C       THIS ROUTINE PERFORMS THE RANDOM WALK FOR ALL PARTICLES
   10 CONTINUE
C get material and particle information
*     U = UINC(1)
*     V = UINC(2)
*     W = UINC(3)
      U = VECT(4)
      V = VECT(5)
      W = VECT(6)
      X = 0.0
      Y = 0.0
      ZZ = 0.0
      BLZNT = 1
      WATE = 1.0
      AGE = 0.0
      NREG = 1
      WTBC = 1.0
C Energy MeV -> eV
*      E = EINC * 1.E6
      E = GEKIN*1.E9
C Material number a la GEANT
*      NMED = NCEL
      NMED = NMAT
      NMEM=1
C reset counter of heavy/charged and gamma bank
      NMEMR = 0
      NMEMG = 0
      INALB=0
      IET=0
      EOLD=E
      UOLD=U
      VOLD=V
      WOLD=W
      OLDWT=WATE
      XOLD=X
      YOLD=Y
      ZOLD=ZZ
      BLZON=BLZNT
      MEDOLD=NMED
      OLDAGE=AGE
      I=1
      CALL GTMED(NMED,IMED)
C get total cross-section
      CALL NSIGTA(E,NMED,TSIG,D,LD(LFP32),LD(LFP33))
C       DETERMINE WHICH ISOTOPE HAS BEEN HIT
      CALL ISOTPE(D,LD,LD(LFP10),D(LFP12),LD(LFP16),LD(LFP26),LD(LFP27),
     +            E,TSIG,IMED,IIN,IIM)
C       THE PARAMETER (IIN) IS THE POINTER FOR ARRAYS DIMENSIONED BY
C       (NNUC) AND THE PARAMETER (IIM) IS THE POINTER FOR ARRAYS
C       DIMENSIONED BY (NMIX)
      LD(LFP42+IMED-1)=LD(LFP42+IMED-1)+1
      INEU = 0
      NNEU = 0
      NHEVY = 0
      NGAMA = 0
      NPSTOR = 0
      CALL COLISN(D,LD,
     + LD(LFP20),LD(LFP21),LD(LFP22),LD(LFP23),LD(LFP24),
     + LD(LFP25),LD(LFP26),LD(LFP27),LD(LFP28),LD(LFP29),LD(LFP30),
     + LD(LFP31),D(LFP34),D(LFP35),LD(LFP41),LD(LFP41+NNUC),
     + LD(LFP42),LD(LFP42+MEDIA),LD(LFP42+2*MEDIA),LD(LFP42+3*MEDIA),
     + LD(LFP42+4*MEDIA),LD(LFP42+5*MEDIA),LD(LFP42+6*MEDIA),
     + LD(LFP42+7*MEDIA),LD(LFP42+8*MEDIA),LD(LFP42+9*MEDIA),
     + LD(LFP42+10*MEDIA),LD(LFP42+11*MEDIA),LD(LFP42+12*MEDIA),
     + LD(LFP42+13*MEDIA),LD(LFP42+14*MEDIA),LD(LFP42+15*MEDIA),
     + LD(LFP42+16*MEDIA),LD(LFP42+17*MEDIA),LD(LFP42+18*MEDIA),
     + LD(LFP42+19*MEDIA),LD(LFP42+20*MEDIA),LD(LFP42+21*MEDIA),
     + LD(LFP42+22*MEDIA),LD(LFP45),LD(LFP46),LD(LFP13),
     + LD(LFP35+NQ*NNUC),D(LFP35+2*NQ*NNUC),IIN,IIM)
      CALL BANKR(D,LD,5)
C -------- fill return arrays with generated particles ---------------
C first heavy/charged particles
   20 NPHETC = 0
      NRECOL = 0
      ERMED  = 0.0
      EETOT = 0.0
C -------- store  neutrons -------------------------------------
      INTCAL = 0
C
      ISTOP = 1
      JPA = LQ(JPART-13)
      AGEMNE = Q(JPA+7)
      NGKINE = 0
      DO 30  N=1,NNEU
         CALL GETPAR(IDNEU,N,IERR)
         IF(IERR.EQ.0) THEN
            NGKINE = NGKINE + 1
            TTKIN  = EP * 1.E-9
            PGEANT = SQRT(TTKIN*(TTKIN+2*AGEMNE))
            GKIN(1,NGKINE) = UP*PGEANT
            GKIN(2,NGKINE) = VP*PGEANT
            GKIN(3,NGKINE) = WP*PGEANT
            GKIN(4,NGKINE) = TTKIN + AGEMNE
            GKIN(5,NGKINE) = 13
            TOFD(NGKINE)   = AGEP * 1.E-9
            GPOS(1,NGKINE) = VECT(1)
            GPOS(2,NGKINE) = VECT(2)
            GPOS(3,NGKINE) = VECT(3)
*           NPHETC = NPHETC + 1
*           IF(NPHETC.GT.MXCP) NPHETC=MXCP
*           IPCAL(NPHETC) = 1
C kinetic energy in MeV
*           EKINET(NPHETC) = EP * 1.E-6
*           UCAL(NPHETC,1) = UP
*           UCAL(NPHETC,2) = VP
*           UCAL(NPHETC,3) = WP
*           CALTIM(NPHETC) = AGEP
         ENDIF
   30 CONTINUE
C -------- store heavy recoil products ------------------------
      DO 40  N=1,NHEVY
         CALL GETPAR(IDHEVY,N,IERR)
         IF(IERR.EQ.0) THEN
C check particle type
            MA = NINT(AMP)
            MZ = NINT(ZMP)
            IF(MA.LE.4.AND.MZ.LE.2) THEN
               IF(NGPART(MA,MZ).EQ.-1) GOTO 40
            ELSE
C get heavy recoil nucleus
               NRECOL = NRECOL + 1
               AMED(NRECOL) = AMP
               ZMED(NRECOL) = ZMP
               ERMED = ERMED + EP * 1.E-9
               GOTO 40
            ENDIF
C store particle type
            NGKINE = NGKINE + 1
            JPA = LQ(JPART-NGPART(MA,MZ))
            AGEMAS = Q(JPA+7)
            TTKIN  = EP * 1.E-9
            PGEANT = SQRT(TTKIN*(TTKIN+2*AGEMAS))
            GKIN(1,NGKINE) = UP*PGEANT
            GKIN(2,NGKINE) = VP*PGEANT
            GKIN(3,NGKINE) = WP*PGEANT
            GKIN(4,NGKINE) = TTKIN + AGEMAS
            GKIN(5,NGKINE) = NGPART(MA,MZ)
            TOFD(NGKINE) = AGEP * 1.E-9
            GPOS(1,NGKINE) = VECT(1)
            GPOS(2,NGKINE) = VECT(2)
            GPOS(3,NGKINE) = VECT(3)
*           NPHETC = NPHETC + 1
*           IF(NPHETC.GT.MXCP) NPHETC=MXCP
*           IPCAL(NPHETC) = NPART(MA,MZ)
C kinetic energy in MeV
*           EKINET(NPHETC) = EP * 1.E-6
*           UCAL(NPHETC,1) = UP
*           UCAL(NPHETC,2) = VP
*           UCAL(NPHETC,3) = WP
*           CALTIM(NPHETC) = AGEP
         ENDIF
   40 CONTINUE
* Number of produced particles (may be > MXGKIN)
      NNEHEG = NGKINE + NGAMA
C
C----------- get generated gammas --------------------
      NS   = 0
      NREM = 0
      DO 50  N=1,NGAMA
         IF (NS.GE.NGAMA) GO TO 60
         NS = NS + 1
         CALL GETPAR(IDGAMA,NS,IERR)
         IF(IERR.EQ.0) THEN
            IF (NNEHEG-NREM.GT.MXGKIN) THEN
               NREM  = NREM + 1
               UP1   = UP
               VP1   = VP
               WP1   = WP
               EP1   = EP
               AGEP1 = AGEP
               NS = NS + 1
*    Get the other gamma to be summed with the previous one
               CALL GETPAR(IDGAMA,NS,IERR)
               IF(IERR.EQ.0) THEN
                 UP  = (UP1*EP1+UP*EP)
                 VP  = (VP1*EP1+VP*EP)
                 WP  = (WP1*EP1+WP*EP)
*    Normalize the new direction cosines
                 WUP = SQRT(UP**2+VP**2+WP**2)
                 UP  = UP/WUP
                 VP  = VP/WUP
                 WP  = WP/WUP
                 EP  = EP1 + EP
                 AGEP = AGEP1 + AGEP
               ENDIF
            ENDIF
            NGKINE = NGKINE + 1
            PGEANT = EP * 1.E-9
            GKIN(1,NGKINE) = UP*PGEANT
            GKIN(2,NGKINE) = VP*PGEANT
            GKIN(3,NGKINE) = WP*PGEANT
            GKIN(4,NGKINE) = PGEANT
            GKIN(5,NGKINE) = 1
            TOFD(NGKINE) = AGEP * 1.E-9
            GPOS(1,NGKINE) = VECT(1)
            GPOS(2,NGKINE) = VECT(2)
            GPOS(3,NGKINE) = VECT(3)
*           NG = NG + 1
*           NPHETC = NPHETC + 1
*           IF(NPHETC.GT.MXCP) NPHETC=MXCP
*           IPCAL(NPHETC) = 11
*           EKINET(NPHETC) = EP*1.E-6
*           UCAL(NPHETC,1) = UP
*           UCAL(NPHETC,2) = VP
*           UCAL(NPHETC,3) = WP
*           CALTIM(NPHETC) = AGEP
C nucleus is in ground state !
            EXMED = 0.0
         ENDIF
   50 CONTINUE
* only one neutron generated -> the particle is the same
   60 IF (NGKINE.EQ.1.AND.GKIN(5,1).EQ.13) THEN
         NGKINE  = 0
         CALL GETPAR(IDNEU,1,IERR)
         VECT(4) = UP
         VECT(5) = VP
         VECT(6) = WP
         GEKIN   = EP * 1.E-9
         GETOT   = GEKIN + AGEMNE
         VECT(7) = SQRT(GEKIN*(GEKIN+2.*AGEMNE))
         TOFG    = TOFG + AGEP * 1.E-9
         ISTOP   = 0
      ENDIF
*
      IF (MTP        .EQ.         2) THEN
         INTCAL = 13
      ELSEIF (MTP    .EQ.        18) THEN
         IF (NHEVY.GT.0) INTCAL = 15
      ELSEIF (MTP    .LT.       100) THEN
         IF (NNEU .GT.0) INTCAL = 20
      ELSEIF (MTP    .EQ.       102) THEN
         IF (NGAMA.GT.0) INTCAL = 18
      ELSEIF (MTP    .GE.       100) THEN
         IF (NHEVY+NGAMA.GT.0) INTCAL = 16
      ENDIF
      IF(NNEU+NHEVY+NGAMA.GT.0.AND.INTCAL.EQ.0) INTCAL=12
      KCASE = NAMEC(INTCAL)
      END
